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Free energies of molecules can be calculated by quantum computations or by normal mode classical 
calculations. However, the first can be computationally impractical for large molecules and the 
second is based on the assumption of harmonic dynamics. We present a novel, accurate and complete 
calculation of molecular free energies in standard classical potentials. In this method we transform 
the molecule by relaxing potential terms which depend on the coordinates of a group of atoms in 
that molecule and calculate the free energy difference associated with the transformation. Then, 
since the transformed molecule can be treated as non interacting systems, the free energy associated 
with these atoms is analytically or numerically calculated. This two-step calculation can be applied 
to calculate free energies of molecules or free energy difference between (possibly large) molecules in 
a general environment. We suggest the potential application of free energy calculation of chemical 
reactions in classical molecular simulations. 
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INTRODUCTION 

Free energy calculations have a variety of applications 
which include binding, solvation and chemical reactions. 
For the applications of solvation and binding equilib¬ 
rium and non-equilibrium methods are used. In equilib¬ 
rium methods (alchemical transformations) a molecule is 
transformed into another through non realistic interme¬ 
diate systems. The average of the derivative of the inter¬ 
polating Hamiltonian is calculated for each intermediate 
system and by numerically integrating, the free energy 
associated with the transformation is calculated. The 
free energies of these transformations are then used to 
calculate the free energy of the desired molecular process 
m- In non equilibrium methods the work in the pro¬ 
cess of switching between the two Hamiltonians is mea¬ 
sured ei n\. For the application of chemical reactions the 
free energy difference between the reactants and prod¬ 
ucts needs to be calculated. Free energies of molecules 
can be calculated by quantum computations |8] or by 
normal mode classical calculations 0 (see Entropies Sec¬ 
tion). While QM calculations are highly accurate they 
are computationally very demanding for large molecules. 
Normal mode calculations can be very efficient but are 
based on harmonic dynamics m- 

Molecular modeling includes covalent bond, bond an¬ 
gle, dihedral angle, improper dihedral angle, electrostatic 
and vdW potentials (see HMSI and Supplementary Ma¬ 
terial). Covalent bond, bond angle and dihedral angle 
terms depend on the coordinates of two, three and four 
nearest covalently linked atoms respectively. Electro¬ 
static and vdW potentials are between every atom pair 
in the system. 

In previous studies, the free energy associated with 
the covalent bond and bond angle terms in molecules 
and dihedral angle term in the context of restraints were 
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directly calculated by assuming that they are harmonic 
(quadratic) and by using the rigid rotor approximation 
and the HJR technique OHS]. However, the free en¬ 
ergies associated with dihedral potential terms in the re¬ 
alistic non-harmonic potentials and with non-harmonic 
covalent bond and bond angle terms m have not been 
calculated. In addition, it has not been suggested to cal¬ 
culate molecular free energies of submolecules with cou¬ 
pled bond angle degrees of freedom such as the methyl 
group. 

Here we present a direct free energy calculation which 
is exact and does not require the potential terms to be 
harmonic (see also Ref. [HI)- We show that the partition 
function associated with a dihedral angle in molecules can 
be decoupled from the partition function of the rest of 
the system independently of the potential function (even 
though the dihedral potential term depends on coordi¬ 
nates of atoms in the system) and calculate analytically 
free energies of realistic dihedral terms. Moreover, we 
show that free energy associated with submolecules with 
coupled bond angle degrees of freedom such as the methyl 
group can be calculated accurately. 

In Section I we present the first step of the calculation 
in which we transform the molecule by relaxing some po¬ 
tential terms of a group of atoms in that molecule and 
calculate the free energy difference associated with the 
transformation. In Section II we present the second step 
of the calculation in which the free energy associated with 
the remaining potentials (which depend on the relative 
spherical coordinates) of this group of atoms is calculated 
analytically or numerically. In Section III we suggest the 
potential application of free energy calculation of chem¬ 
ical reactions. The demonstartions of the method are 
explained in the Demonstration section. 

Alchemical free energy transformation (as in Section I) 
is a standard procedure in molecular free energy calcu¬ 
lations mmmm and is usually used in the context 
of solvation and binding or in the context of chemical re¬ 
actions in combination with QM calculations )8j. Here 
we apply it in the context of free energy calculations 
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of molecules or submolecules in a general environment, 
where the second calculation is used to obtain the free en¬ 
ergy difference between two molecules with submolecule 
in common. We combine the calculation with novel direct 
free energy calculations explained in Section II, followed 
by the novel application of free energy of chemical reac¬ 
tions presented in Section III. Here we do not need to 
assume harmonic dynamics. 

I. THE TRANSFORMATION 

We now explain the first step of the calculation in which 
we we transform the molecule by relaxing potential terms 
that depend on the coordinates of a group of atoms in 
that molecule and calculate the free energy difference as¬ 
sociated with the transformation using Thermodynamic 
Integration (TI) [3][2TJ[22]. Alternatively, it can be per¬ 
formed using methods such as Exponential Averaging/ 
Free Energy Perturbation (FEP) [23], Bennett Accep¬ 
tance Ratio (BAR) and Weighted Histogram Anal¬ 
ysis Method [25] (WHAM). We denote the Hamiltonian 
with the terms that are removed in the transformation 
by H r and the Hamiltonian with the other terms by H c . 
We define the Hamiltonian H r as the improper dihedral, 
electrostatic and vdW terms that depend on the coordi¬ 
nates of the group of selected atoms. The A dependent 
Hamiltonian that transforms between the systems can be 
written as follows: 

H (A) = XH r + H c , 

where A G [0,1]. The free energy difference between the 
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FIG. 1. Systems in the transformation of phenol at A — 
0,0.5,1. The semi transparent/ transparent atoms represent 
atoms with interactions partly/ totally removed 

original and transformed systems can be written as: 

AF a ^ a , = ~k B T^ rfln ^ A (/? ’ A) riA = £ (H r ) dX, 

where A' is the transformed system in which A = 0. Sim¬ 
ulating intermediate systems at A values in the range [0,1] 
and calculating (JT r ), followed by numerical integration, 
will enable us to calculate the free energy difference. It is 
noted that this calculation is valid for any environment 


(gas, solvent etc.) since the interactions of the group 
of atoms with the environment are also relaxed in the 
transformation. In addition, soft core potentials are used 
in the Vdw and electrostatic interactions to avoid diver¬ 
gences of (H r ) at A —>> 0 [26] . 

As an example we present in Fig. [ljthree systems in the 
transformation of the molecule phenol with the selected 
atoms O and HI. The Hamiltonian H r includes the vdW 
and electrostatic interactions of the atoms O and HI with 
all the atoms in the system, and the improper dihedral 
term 0 ci-C 2 -o-C 3 * 


II. THE REMAINING FREE ENERGY 
DIFFERENCE 

We now turn to explain how the free energy associated 
with the selected atoms can be calculated. We first switch 
to relative coordinates and then to spherical coordinates. 

n k n 

dft = II d r> i = rfr'i dri 

i= 1 i=2 i=k-\-1 

k n 

= dr ' 1 n dri n r i sin 

i—2 i=k-\-l 

where = r' — r'_ x , which will be chosen as the position 
of atoms relative to a covalently bound atom, k + 1 rep¬ 
resents the first atom in the group of selected atoms (e.g 
the atom O in the example of phenol). 

The transformed molecule A’ is first divided into ele¬ 
ments of standard covalent bonds, bond junctions, dihe¬ 
dral terms and of more complex structures that include 
molecular rings. The dihedral potential term depends 
on the angle between two planes which is equal to the 
< j) angle in spherical coordinates. This correspondence 
can be understood by recalling the definition of the an¬ 
gle between two planes as the angle between the vectors 
in these planes that are perpendicular to the intersection 
line of these planes. Hence, the covalent bond, bond an¬ 
gle and dihedral angle terms can be expressed in terms 
of the spherical variables r, 6 and 0 defined with respect 
to the relevant atoms. 

For the example of phenol we can write the partition 
function of transformed phenol as follows 

Za> = f e ~ pH ^ ) e~ l3Hr o( r o) e -0 H r- H (rH) x 

e -pHe H (o H ) s i n Q He -P H 4> H (<t>H)d^rQdrorjjdrHdOndfpH, 

where the degrees of freedom of atom H and O are de¬ 
noted by rn,#H,0H and ro,#o,0o respectively. All the 
degrees of freedom of all the other atoms in the system 
(molecule+solvent molecules in the case of explicit sol¬ 
vent) and #0,00 are denoted by ft. We define Oh and 
4>h as follows: Oh = 0hi-o-C2,4>h = 4>hi-o-C2~ci- 
Our goal in this example is to calculate the free energy 
associated with all the degrees of freedom except for the 
ones denoted by ft. The free energy of the subsystem 



3 


with the degrees of freedom denoted by ft can for exam¬ 
ple then cancel out with an identical subsystem when cal¬ 
culating the free energy difference between two molecules 
with a submolecule in common. Note that for example 
4>h does depend on the coordinates of the Cl atom which 
is included in ft and therfore the decomposition of the 
partition functions is not trivial. 

As will become clear, the integration in each element 
is independent of the others. Thus the integrals can be 
performed separately and then multiplied to yield the 


partition function and hence the free energy difference. 

We write the free energy factors explicitly: 

Covalent bond 

The partition function of a covalent bond represented 
by a harmonic (quadratic) potential can be written as 
follows (for non-quadratic terms integrated analytically 
see Ref. El): 



/ 


f3k c (r—d) 2 2 

e 2 r dr 


(2d 2 f3k c + l) (erf (dy//3k c ) + l) + 2^/ire d kcdj2 d^/ [3k c 

p (/?fcc ) 3/2 


( 2 ) 


where l is an arbitrary length ( l 3 cancels out in compar¬ 
isons). To account for the bond-dissociation energy/bond 
energy one has to differentiate between the potential at 
r = 0 and k c 0 (free atoms) either by introducing a 
factor 27 29; or by using potentials such as the Morse 
potential (numerical integration) [3]3j. 


the rest of the molecule moves as a rigid body. Since the 
spherical coordinates representation includes three inde¬ 
pendent variables, varying a bond angle is decoupled from 
all the other degrees of freedom of the molecule. Hence 
the calculation of free energy factor associated with the 
bond angle potential term is straightforward. For a har¬ 
monic (quadratic) term we can write : 


Two Bonds Junctions 

When considering the case of a Linear/Bent molecular 
shapes, it can be seen that when varying the bond angle, 



Z b = J e-M sin OdO = / e-^ ^ 
9 0 /3ke + (dkeir 
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i Odd = 
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e ie ° X 


ierfi 


1 + i/3ko( 7T - Op) 

7We 


(3) 


This expression is real for positive and real values of ko : (3 
and Oq. 

When varying a dihedral angle, the potential term 
value depends on the orientation of first bond (which de¬ 
termines the axis from which the dihedral angle is mea¬ 
sured). However, since the integration has to be per¬ 
formed over all the range [0, 27r], integrating with respect 
to the 0 angle will yield a factor which is independent 
of the location of the first bond. Thus, the integration 
does not depend on the direction of the first bond and is 
straightforward. For the commonly used dihedral angle 
potential we get 

Zd = j e-^P+^^-^dcp = 2ire~ l3k *I 0 (/?fc 0 ) ,(4) 

where Iq is the modified Bessel function of the first kind. 
Note that this result does not depend on n, which is an 
integer. In the context of restraints quadratic dihedral 
terms can be used and we give the exact expression for 


the corresponding partition function in the Appendix for 
completeness. 

It can now be seen that the partition function of trans¬ 
formed phenol in Eq. (1) can be decomposed (exactly) 
into the partition function associated with the ft degrees 
of freedom and the covalent bond, bond angle and dihe¬ 
dral angle partition functions. 


Three or more Bonds Junctions 

Molecule shapes can include a monomer (covalent 
bond) that splits into more than one monomer. Such ex¬ 
amples are the trigonal planar, tetrahedral trigonal pyra¬ 
midal etc. These cases will necessitate numerical integra¬ 
tion which can be performed using the Spherical law of 
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cosines: 


cos (# 12 ) = cos ( 0 i) cos ( 02 ) + sin ( 0 i) sin ( 62 ) cos (A 0 ), 

( 5 ) 

where 0 i, 62 denote the bond angles of two bonds with 
respect to the principal bond and 0 i 2 , A 0 denote the bond 
angle and the difference in 0 angle between these two 
bonds respectively. Usually in these cases there is one 
dihedral term, which we denote by 0 1 , that depends on 


the angle defined by one of the bonds, the principal bond 
and a previous bond. Since the integration over the other 
degrees of freedom yields a factor that is independent of 
the value of 0i, the integrations are decoupled. Thus, the 
integration for the case of one monomer that splits into 
two can be written as follows: 

Z = J e -«v»+v , «Od0idM0 1 #2 = Z d Z 3b , (6) 

where Z& is according to the definition in Eq. @ and 


Z 3b = Je 9 °f+ k ^ 0 2) 2 +fci2(0i2 s m(0 2 )d9 1 d9 2 dth. 

For the general case it can be written as follows: 


/ n n 

© smOidOi n #», 

i i>j i i> 2 


( 7 ) 

(8) 


where 0^- can be calculated from Eq. © and 0 i, which 
has to be substituted in A 0 in this equation, is an ar¬ 
bitrary constant. In case there are energy terms that 
introduce complexity they can be relaxed in the trans¬ 
formation. For an explicit example with a bond junction 
see SM. 

The free energy factors can be substituted in: 

= (X> Z ,„ +2>Z*, + £.n Z J , 

V i i i / 

(9) 

to give the free energy difference, where A" denotes a 
submolecule of the transformed molecule that is not nec¬ 
essarily realistic whose free energy will not be calculated. 

In addition, the partition function of complex struc¬ 
tures may be calculated numerically. In many cases the 
complex structures need to be compared to identical ones, 
eliminating the need for these calculations. 

Thus, we can write in terms of the partition functions: 

Imp q r 

Z —>• Z[ nt n z n z ' n %2bi ^3 bi n ^complex^ •> 

i=1 i=l i= 1 i=l i= 1 

where Z[ nt represents the partition function of the sub¬ 
molecule that is fully interacting and Z Ci and Z^. repre¬ 
sent the ith covalent bond and dihedral angle partition 
function respectively. Z 26 . and Z 35 . represent the ith 
two bond and three or more bond junctions respectively 
and Complex, represents the ith complex structure parti¬ 
tion function. The arrow represents the transformation 
A = 1 0 . 

We can also choose the group of selected atoms as all 
the atoms of the molecule, and calculate the free energy 
in a similar manner. The partition functions in this case 
can be written as follows: 

Imp q r 

Z Zfp n^n Zdi %2bi Z 3 bi n ZcompX ex^ 5 

i=1 i=l i=1 i=1 i=1 


where Zf p denotes the partition function of a free particle. 

In Fig. [2] we present the free energy decomposition of 
transformed phenol (A’). The degrees of freedom asso¬ 
ciated with each system are written on top. F denotes 
the free energy and A” which is the reference (possibly 
imaginary) submolecule, is presented in the first term 
in the decomposition. It should be noted that the 0o 
and 0 o degrees of freedom are associated with the sub¬ 
molecule A” and the ro degree of freedom is associated 
with the second system in the decomposition. The po¬ 
tential terms that depend on the coordinates of the O 
and H atoms in the five decomposed systems are: 1 : 
Vb (0Cl-C2-o) J Vb (@C3-C2-o) ? Vd (0C4-C1-C2-O) 2 : 
V c (ro) , 3 : V c (rn ), 4 : Vb(0c2-o-Hi) and 5 : 

(0ci-C2-o-izi)- The reason that the degrees of free¬ 
dom 0q and 0 o are associated with the first system in the 
decomposition is that the corresponding potential terms 
usually do not depend on the atom (e.g O). The free 
energies of the second and third systems are calculated 
according to Eq. © and the free energies of the fourth 
and fifth systems are calculated according to Eqs. © 
and © respectively. 


DEMONSTRATIONS 

The covalent bond free energy calculation has been 
demonstrated and compared with numerical integration 
for two molecules of two atoms in a spherical poten¬ 
tial (see Supplementary Material). In addition, for 
methanethiol molecule the free energies were calculated 
according to Eqs. (2)-(7) and were in agreement with 
the results obtained by performing transformations in 
MD simulations. The computations using these equa¬ 
tions were faster by factors varying between 5 • 10 3 and 
10 12 (for the same computational resources) p~7 . 
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FIG. 2. The decomposition of the free energies of the trans¬ 
formed molecule A’ 


III. POTENTIAL APPLICATION 

Here we suggest the potential application of calculating 
free energies of chemical reactions using classical molec¬ 
ular simulations followed by analytic or numerical calcu¬ 
lations. We will calculate the free energies of only small 
parts of the molecules, possibly in solvent, to obtain the 
free energy of the chemical reaction which can involve 
large molecules. 

To obtain the equilibrium constant of a chemical reac¬ 


= AFA—tA' + - AFj 


An example for such a calculation is given in Fig. [3] 
The molecules that participate in the chemical reaction 
and the transformed molecules are presented on top and 
bottom respectively. In this case molecule B is matched 
with C and Fb>> = FcF a" and Fd", which are the 
free energy of a free particle, are equal. 

We presented an accurate and complete method for 
calculating molecular free energies in classical potentials. 
We suggested the potential application of free energy cal¬ 
culation of chemical reactions in classical potentials. This 
application can find use in organic chemistry, biochem¬ 
istry and drug discovery. 
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tion the standard Gibbs free energies (which are similar 
to Helmholtz free energies in many cases m) are usu¬ 
ally calculated. The standard state is the hypothetical 
state with the standard state concentration but exhibit¬ 
ing infinite-dilution behavior (the interactions between 
e.g the solute molecules are negligible). Hence, when we 
are interested in the standard free energy of one substance 
a single copy of the molecule of interest can be simulated 
either in vacuum or solvent environments. 

The idea is to transform the reactants and the prod¬ 
ucts, between which the free energy difference is calcu¬ 
lated, into molecules that have the same partition func¬ 
tions up to factors that can be calculated. First, we 
match reactant molecules with product molecules that 
have submolecules in common if possible so that the free 
energy of these submolecules will cancel out. Then, we 
transform the molecules by relaxing potential terms of 
the atoms that are different between the molecules and 
calculate the free energy differences associated with each 
transformation. Finally the free energy factors associ¬ 
ated with the different atoms are calculated and we can 
deduce the free energy of the chemical reaction. For ex¬ 
ample, given the chemical reaction 

A + B^C + D, ( 10 ) 

we can match, if possible, molecules A , B to molecules 
C, D. Then we transform each of the molecules to 
A'^B'^C' and D' and calculate the free energy dif¬ 
ferences AF a ^A',AF b ^b',AF c ^c' and A F D ^ D > re¬ 
spectively. Finally we calculate the free energy factors 
AFa'^a" , AFb'^b" , AFc'^c" and A F D >^ D n as previ¬ 
ously explained. Thus we get 


»' + AFa'^a" + A F B '^b" ~ AF C '^c" ~ AF D >^ D n. (11) 


Appendix 

The exact expression for the partition function associ¬ 
ated with a quadratic dihedral term is the following 

^quadratic = f = 


VI [erf ((ftpygfcd) ~ erf ((</> 0 - 2 tt) y/ffid)\ 

2 VWd 
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FIG. 3. A scheme of the free energy calculation of the 
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